Estimating Granger causality from Fourier and wavelet transforms of time series data 
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Experiments in many fields of science and engineering yield data in the form of time series. 
) , The Fourier and wavelet transform-based nonparametric methods are used widely to study the 

^— ( ' spectral characteristics of these time series data. Here, we extend the framework of nonparametric 

' spectral methods to include the estimation of Granger causality spectra for assessing directional 

influences. We illustrate the utility of the proposed methods using synthetic data from network 
models consisting of interacting dynamical systems. [Physical Review Letters, in press]. 
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Extracting information flow in networks of coupled dynamical systems from the time series measurements of their 
activity is of great interest in physical, biological and social sciences. Such knowledge holds the key to the understand- 
ing of phenomena ranging from turbulent fluids to interacting genes and proteins to networks of neural ensembles. 
Granger causality has emerged in recent years as a leading statistical technique for accomplishing this goal. The 
definition of Granger causality [l[ is based on the theory of linear prediction 0] and its original estimation framework 
requires autoregressive (AR) modeling of time series data fH ^ . Such parametric Granger causality and associated 
spectral decompositions have been applied in a wide variety of fields including condensed matter physics 4] , neuro- 
S i H genetics i, climate science [13, El, and economics |12|. However, the parametric modeling 



CIh, science 

methods often encounter difficulties such as uncertainty in model parameters and inability to fit data with complex 
spectral contents [l3| . On the other hand, the Fourier and wavelet transform-based nonparametric spectral methods 
^ ■ are known to be free from such difficulties and have been used extensively in the analysis of univariate and 
0^ multivariate experimental time series A weakness of the current nonparametric framework is that it lacks 

the ability for estimating Granger causality. In this Letter, we overcome this weakness by proposing a nonparametric 
, approach to estimate Granger causality directly from Fourier and wavelet transforms of data, eliminating the need 
• ' of explicit AR modeling. Time-domain Granger causality can be obtained by integrating the corresponding spectral 
representation over frequency [3]. Below, we present the theory and apply it to simulated time series. 

Grar^ger causalUy: the parameMc estirr^at^or. approach. Granger causality Q is a measure of causal or directional 
, influence from one time series to another and is based on linear predictions of time series. Consider two simultane- 
IL^ ously recorded time series: Xi : a:i(l), a;i(2), ...,xi{t), X2 : 2:2(1), X2(2), ...,X2{t), ... from two stationary stochastic 
processes {Xi,X2). Now, using AR representations, we construct bivariate linear prediction models for xi{t) and 
X2it): 



xiit) = ^bii,]Xi{t - j) +^bi2,jX2{t - j) + eii2{t) (1) 

00 00 

X2{t) = X!^2ija;i(i- j) +^622jX2(i-j) +e2|i(i) (2) 

along with the univariate models: xi{t) = X^jli <^jXi{t — j) + ei{t) and X2{t) = J^'jLi PjX2it — j) + e2(i)- Here, e's 
are the prediction errors. If var(ei|2(t)) < var(ei(t)) in some suitable statistical sense, then X2 is said to have a causal 
influence on Xi. Similarly, if var(e2[i(i)) < var(e2(t)), then there is a causal influence from Xi to X2. These causal 

influences are quantified in time domain 31 by Fj^i = In 7^—^77, where i = 1, 2 and j — 2, 1. 

\ ' var(e,|j(t)) 

Experimental processes are often rich in oscillatory content, lending themselves naturally to spectral analysis. 
The spectral decomposition of Granger's time-domain causality was proposed by Geweke in 1982 [3|]. To derive the 
frequency-domain Granger causality, we start with Eq. ((T][2]) . We rewrite these equations in a matrix form with a lag 
operator L: Lx{t) = x{t — 1) as 



bn{L) bi2{L) \ { xi{t)\ ^ /ei|2 

b2l{L) h22{L) j I X2{t) j U2I1 



(3) 
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where bij{L) = X]fc°=o ^ij.kL^ with bij^ = Sij (the Kronecker delta function). The covariance matrix of the noise terms 
is S = ( ) where Sn = var(ei|2), S12 = ^21 = cov(ei|2, €211), and S22 — var(e2|i)- Fom'ier transforming 



^21 ^22 
Eq. (131) yields 



(4) 



where the components of the coefficient matrix [Bij{f)] are Bimif) = Sim — J2'k'=i bim.,ke jj-^ terms of transfer 

fmiction matrix (H(/) = [Bij{f)]^^), Eq. (g]) becomes 



X2{f) \H2l{f) H22U) \E2if) 



(5) 



Then, the spectral density matrix S(/) is given by 

S(/) - H(/)SH*(/), (6) 

where * denotes matrix adjoint. To examine the causal influence from X2 to Xi, one needs to look at the auto-spectrum 
of xi(t)-series, which is 6*11 (/) = HiiT^nHli + 2Ei2Re(_ffii_ffj*2) + i7i2S22-ffr2- Here, because of the cross-terms in 
this expression for S'n, the causal power contribution is not obvious. Geweke Q introduced a transformation that 
eliminates the cross terms and makes an intrinsic power term and a causal power term identifiable. For Xi-process, 

this transformation is achieved by left-multiplying Eq. (j4]) on both sides with ( v ^/v ? )' "^hich yields: 



-S12/S11 1 



i?ii(/) B,2U) \ ( X^{f) \ ^( E,U) 

B2l{f) B22{f) J \ X2{f)) \E2{f) 



(7) 



where B^iiJ) = B^iif) - B22U) = B22U) - |i^Bi2(/), and E2U) = E2{f) - ^Ei{f). The elements 

Zjii Zjli Zjli 

of the new transfer function H(/) then become 7?n(/) = iJn(/) + ^H^^if), HMf) = Huif), H21U) = H21U) + 

^11 

S12 ~ ~ ~ Y? 

— — Hii{f), and H22{f) — H22if)- Here, cov(Ei, E2) = and the new variance of X2{t) is S22 = S22 — ■=r^. Now, the 

auto-spectrum of xi{t) is decomposed into two obvious parts: Sii{f) — + Hi2{f)T,22Hi2{f)j where 

the first term accounts for the intrinsic power of xi{t) and the second term for causal power due to the influence from 
X2 to Xi . Since Granger causality is the natural logarithm of the ratio of total power to intrinsic power [3] , causality 
from X2 to Xi (or, 2 to 1) at frequency / is 

/2^i(/) = In . ^"^{.l X ' (8) 

Siiif) - (Y22 - \Hi2iW 

using the expressions for Sn and E22 obtained after the transformation. Next, by taking the transformation matrix 
1 — S12/S22 
1 
ission for 

time-domain measure is theoretically related to the frequency-domain measure as i^2^i < — / l2^i{f)df , but for 

all processes of practical interest, the equality holds. 

From the above discussion, it is clear that the estimation of frequency-domain Granger causality requires noise 
covariance and transfer function which are obtained as part of the AR data modeling. The mathematics behind 
this parametric approach to obtain these quantities is well-established. However, for nonparametric methods the 
current estimation framework does not contain provisions for computing these quantities. Moreover, the parametric 
estimation method from finite data can often produce erroneous results if the series in Eq. ((l][2|) are not truncated 
to proper model orders. There are criteria [l^ for choosing proper AR model order, but these criteria cannot always 
be satisfied. In addition, AR modeling approach does not always capture all the spectral features To overcome 
these difficulties, we propose a nonparametric estimation approach, in which we derive, based on Fourier and wavelet 



and performing the same analysis, one can get Granger causality /i_>2(/) from Xi to X2, the 

expression for which can be obtained just by exchanging subscripts 1 and 2 in Eq. Geweke[3| showed that the 
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transforms of time series data, noise covariance and transfer function to be used in Geweke's formulae such as Eq. ([5]) 
for Granger causality estimates. 

Granger causality: the nonparametric estimation approach. In the nonparametric approach, spectral density ma- 
trices are estimated by using direct Fourier and wavelet transforms of time series data. These matrices then un- 
dergo spectral density matrix factorization [l^i 12 1 and Geweke's variance decomposition [5] for the estimation of 



Granger causality. To explain this approach, let us consider a bivariate process with time series xi{t) and X2{t), 
their Fourier transforms: Xi{f) and X2{f), and wavelet transforms: Wxi{t,f) and Wx2{t,f)- Then, the spectral 

? S S \ 

matrix S is defined as: S = „ „ ), where, in the Fourier-based method, Sim — {Xi{f)X„i{f)*) , and, in 

\ 'J21 '-'22 / 

the wavelet method, Sim ~ {^Xi{t, f)Wx„^{t, f)*) ■ Here, I — I, 2, m — I, 2, and (.) is averaging over multiple 
realizations. Smoother Fourier-based spectral density with reduced estimation bias can be obtained by using the 
multitaper technique which involves the use of discrete spheroidal sequences (DPSS) [2l[. The continuous 

wavelet transform Wx, (t, s) at time t and scale s is computed by the convolution of time series xi with a scaled and 
translated version of a prototype wavelet ^'(77) that satisfies zero-mean and unity square- norm conditions (23 , H^: 

Wx, (t, s) = —= / drj'i/* I xi{rj). Using the relationship between s and / for a given prototype wavelet, such 

V^i-oo As/ 

as the Morlet wavelet [23|, |2J|, one can transform Wxi (t, s) into Wxi {t, /). The wavelet transform at / = can be ob- 
tained by a numerical extrapolation. S(/) or S(t, /) thus formed is a square matrix that can be defined in the interval 
[—IT, tt] and, for all processes of practical interest, satifies the following properties: (i) S{d) is Hermitian, nonnegative, 
and S{—6) — 5^(6*), where 6 = 27r/ and ^ denotes matrix transpose, and (it) S{9) is integrable and has a Fourier series 
expansion: S{9) = J2^=-oo Ike^''^ , where the covariance sequence {jk}'^oo is formed by jk = (l/27r) J^^ S{9)e~'^^^d9. 

According to the factorization theorem [I^, 25 1, the spectral density matrix S as defined above can be factored [2^ 
into a set of unique minimum-phase functions: 

S = VV'*, (9) 

where * denotes matrix adjoint, ip{e^^) — J2T=q ■^'^^^''^ defined on the unit circle {\z\ — 1}, and = 
{l/2Tr) J^^iP{e''')e-'''''d9. Moreover, V 

can be holomorphically extended [2^ to the inner disk {|z| < 1} as 

ip{z) — '^'kLo ^kz'^ where tp{0) = Aq, a real, upper triangular matrix with positive diagonal elements. Similarly 
S and H can be defined as functions of z with H{0) — I. Comparing the right hand sides of Eqs jS]) and ^ at z = 
we get 

S = AoAj. (10) 

Rewriting Eq. ([5]) as S = AqA'^Aq "^V^* and comparing with Eqs (O and pHl) . we arrive at the expression for 

the transfer function: 

H = Mo' (11) 

Now, by substituting the specific elements of the noise covariance and transfer function from Eq. (jlOp and (|lip into 
Eq. ([5]), one can estimate pairwise Granger causality spectra. In case of the wavelet, these calculations are repeated 
along the time axis for each time point to get the time-frequency representation of Granger causality. 

Numerical examples. We consider network models with two autoregressive processes Xi and X2 as nodes where 
Xi{t) = 0.55Xi{t - 1) - 0.8Xi{t - 2) + C{t)X2{t - 1) + e{t) and X2(t) = 0.55X2(i - 1) - 0.8X2(t - 2) -h ^(t). Here, t 
is a time index, e{t) and S,{t) are independent white noise processes with zero means and unit variances, C{t) is the 
coupling strength, and the sampling rate is considered to be 200 Hz. By construction, only X2 has a causal infiuence 
on Xi. First, we fix C{t) at 0.2 V t, generate dataset of 5000 realizations (trials), each consisting of 5000 data points, 
and apply the Fourier-based nonparametric method. The power spectra of Xi and X2 (Fig. 1(a)) and coherence 
spectra between Xi and X2 (Fig. 1 (b)) show 40 Hz peaks. Figure 1 (c) shows the Granger causality spectra. 
Here, both the nonparametric (NP) and parametric (P) approaches yield identical results, recovering correctly the 
underlying directional influences. Since the proper model order is chosen here and the dataset is large, the parametric 
causality estimates can be assumed to represent the theoretical values. Next, we let the unidirectional coupling of 
X2 to Xi change in its strength C{t) over time as shown in Fig 2(a), generate 5000 realizations, each with 900 
time-points. Then, the wavelet spectra are computed for all the trials using the Morlet wavelet (as used in [1^). The 
average wavelet spectra are obtained by averaging over these individual spectra. The average spectra at a time point 
is subsequently factored, and H and S are obtained and used in Eq. ([8]) to obtain Granger causality spectra. By 
repeating these calculations along time axis, one gets the complete time-frequency maps of Granger causality (Fig. 
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FIG. 1: (a) Power, (b) coherence, and (c) Granger causality spectra from both Fourier transform-based nonparametric (NP) 
and parametric (P) methods. There is an exceUent agreement between NP and P estimates. 



2(b-c)), which also recovers the correct directional influences. Granger causality magnitude increases with coupling 
strength. 

Here, the proposed Granger causality techniques are tested on datasets with a large number of long trials. However, 
these techniques can also be used reliably with fewer trials. Increasing the number of trials leads to spectral estimates 
with smaller variance. A single, sufficiently long stationary time series can be broken into smaller segments, each of 
which can be treated as a distinct trial. The use of multitaper and multiwavelet f27j techniques can yield better 
estimates of Granger causality in case of a dataset with shorter length and fewer trials. See Supplementary Material 
[2^ for additional numerical examples and applications to brain signals. 

Conclusion. Granger causality is a key technique for assessing causal relations and information flow among simul- 
taneous time series. Its traditional parametric estimation framework often suffers from uncertainty in model order 
selection and inability to fully account for complex spectral features. We develop a nonparametric approach based on 
the direct Fourier and wavelet transforms of data that eliminates the need of parametric data modeling and extends 
the capability of Fourier and wavelet-based suites of nonparametric spectral tools. It is expected that the integration 
of the proposed method into existing laboratory analysis routines will provide the basis for gain ing deeper insights 
into the organization of dynamical networks arising in many fields of science and engineering [29|. 

We thank G. T. Wilson for useful email communications. This work was supported by NIH Grant MH71620. GR 
was also supported by grants from DRDO, DST (SR/S4/MS:419/07) and UGC (under SAP-Phase IV). 
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FIG. 2: (color online). Wavelet-based Granger causality: time-frequency representation of causality. Fig (a); temporal 
structures of couplings constructed in the network model: the coupling of X2 with Xi stays 0.25 during time [0, 2] sec, slowly 
changes to during [2, 2.25] sec, and stays during time > 2.25 sec, whereas the coupling of Xi with X2 is throughout. The 
slow transition in the middle is modeled by the tangent of a hyperbolic function. Fig (b, c): time-frequency maps of Granger 
causality correctly represent the temporal structures of couplings as in Fig (a) for the network model. 
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